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Abstract 

With the completion of the Planck mission, in order to continue to gather cosmological information it 
has become crucial to understand the Large Scale Structures (LSS) of the universe to percent accuracy. 
The Effective Field Theory of LSS (EFTofLSS) is a novel theoretical framework that aims to develop an 
analytic understanding of LSS at long distances, where inhomogeneities are small. We further develop the 
description of biased tracers in the EFTofLSS to account for the effect of baryonic physics and primordial non- 
Gaussianities, finding that new bias coefficients are required. Then, restricting to dark matter with Gaussian 
initial conditions, we describe the prediction of the EFTofLSS for the one-loop halo-halo and halo-matter 
two-point functions, and for the tree-level halo-halo-halo, matter-halo-halo and matter-matter-halo three- 
point functions. Several new bias coefficients are needed in the EFTofLSS, even though their contribution at 
a given order can be degenerate and the same parameters contribute to multiple observables. We develop a 
method to reduce the number of biases to an irreducible basis, and find that, at the order at which we work, 
seven bias parameters are enough to describe this extremely rich set of statistics. We then compare with the 
output of an A-body simulation where the normalization parameter of the linear power spectrum is set to 
(Tg = 0.9. For the lowest mass bin, we find percent level agreement up to A: ~ 0.3/iMpc“^ for the one-loop 
two-point functions, and up to /c ~ 0.15/i Mpc“^ for the tree-level three-point functions, with the /c-reach 
decreasing with higher mass bins. This is consistent with the theoretical estimates, and suggests that the 
cosmological information in LSS amenable to analytical control is much more than previously believed. 
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1 Introduction 

With the completion of the Planck mission, most of the future cosmological information will probably 
need to come from Large Scale Structures (LSS). It is therefore essential to understand their behavior. 
In particular, it is of great interest to explore until what high fc-value we can push our theoretical 
control. This is important because it gives an indication of how much information we have hope to 
extract about the primordial universe from LSS surveys. 

The Effective Field Theory of Large Scale Structures (EFTofLSS) [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 
13, 14, 15, 16] is a theoretical framework that aims at developing a rigorous analytic understanding 
of the clustering of LSS, through a perturbative expansion in the size of the fluctuations. It does 
so by correctly parametrizing the effect that short distance non-linearities have at large distances 
through an effective stress tensor where some parameters, such as the speed of sound, the viscosity, 
etc, need to be measured either from observations or from numerical simulations. 
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So far, the EFTofLSS has been compared to data only for the case of dark matter correlation 
functions, specihcally the 2- and 3-point functions at redshift ^ = 0 [2, 3, 12, 10, 11], and the 2-pt 
functions has just been studied at all redshift in [16] . In these cases, the results have been extremely 
promising. For the power spectrum, former analytical techniques start failing at fc ~ 0.1hMpc“^ 
while the FFT, with the inclusion of only one single parameter, fails at fc ~ 0.3hMpc“^ at one- 
loop order, and at fc ~ 0.6hMpc“^ at two-loops. Consistently, the dark matter bispectrum at 
one-loop order agrees until k ~ 0.3hMpc“^ by using the same single parameter [10]. Similarly, and 
consistently, the momentum power spectrum has been shown to agree with data up to fc ~ 0.3h Mpc“^ 
at one-loop order, in practice using just the same parameter [9]. The slope of the power spectrum of 
the vorticous component of the velocity has also been correctly predicted [3, 8]. A recent extension 
of the FFTofLSS to two species has shown that baryonic effects can be consistently included in the 
EFT, and this leads to the same UV reach when star formation physics is present as for dark-matter 
only simulations with the addition of just one parameter, the speed of sound of baryons [14]. 

These results are extremely promising because the available number of modes scales like the cube 
of the maximum wavenumber, so that an improvement of a factor of six in the /c-reach translates in 
a factor of about 200 in the number of available modes. Such a gain, if conhrmed for all observables 
in LSS, would revolutionize our expectations for what we are going to learn from LSS about the 
primordial universe. 

Comparably important is the fact that in the EFTofLSS it is possible to estimate the size of the 
terms that are not included in a calculation at a given order, so that one can forecast the UV reach 
of a given calculation, within order one numbers. The results obtained so far with the EFTofLSS 
meet this requirement. Furthermore, since the EFTofLSS represents the description of a dynamical 
system, the same parameters affect many observables, so that the predictive power of the theory is 
still very large even in the event that one were to measure the parameters of the EFT directly from 
observations rather than from simulations. These are characteristics that distinguish a theory, such 
as the EFTofLSS, from a model. 

The purpose of this work is to compare the predictions of the EFTofLSS for halos to numerical 
simulations. The relevant formulas for biased tracers in the EFTofLSS were worked out in [12], 
importantly generalizing and completing the treatment of [20]. The construction of [12] writes down 
all the bias coefficients that connect the distribution of galaxies or halos to the one of dark matter. In 
particular in [12] it was pointed out that the distribution of collapsed objects at a given time does not 
depend on the underlying dark matter held at the same instant of time, but on the whole history of 
the dark matter held. This is forced by the fact that the formation time for a gravitationally collapsed 
object is of order of the Hubble scale, which is the same time scale over which long wavelength modes 
evolve. This has the consequence that the collapsed objects depend on the underlying dark matter 
held through a diherent coefficient according to the order at which the dark matter held is evaluated. 
Equivalently, the collapsed objects depend also on the convective time derivatives of the dark matter 

^By former analytic techniques, we mean Standard Perturbation Theory (SPT). In the literature there exist several 
versions of perturbations theory where a subset of the SPT diagrams is resummed, such as Renormalized Perturbation 
Theory (RPT) [17], or the Lagrangian formulation (see for example [18, 19]). As explained in [12], when correctly 
implemented, these resummation techniques should give the same UV reach as SPT for the considered quantities. 
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field, with the peculiarity that this time derivatives are unsuppressed (as the theory is non-local in 
time). 

As already pointed out in [12], at a given order in perturbation theory, the contribution of some 
of these terms are potentially degenerate. In Sec. 2 we will develop a method to reduce the basis 
of terms in [12] to a set of linearly independent ones (see also [15] where an analogous irreducible 
basis has recently been constructed from the basis in [12] in a way that is different, and maybe more 
elegant, from the one we do here). This is not quite the full simplihcation we will hnd. In fact, some 
terms become degenerate only in some observables. For example, we will hnd that several cubic 
bias terms, which are independent at the level of helds, contribute in a degenerate way to the power 
spectrum. 

In this paper, we will mainly limit ourselves to study at redshift z = 0 the halo-halo and halo- 
dark matter 2-point functions, Phh and Phm, at one-loop order, and all the bispectra constructed 
from the halo held and the dark matter held, Bhhh, B^mm, at tree level. After this study, we 

will hnd that seven bias parameters are necessary to describe these quantities at the order we wish. 
We therefore expect that the UV reach will be about k ~ 0.3hMpc“^ for one loop quantities, and 
k ~ 0.1h.Mpc“^ for tree level quantities. At the order we work, we hnd that the terms we need to 
include are the same that were already present in [20]. This is just an accident of perturbation theory. 
The new terms that are identihed in [12] represent a relevant physical ehect that is not present in the 
formalism of [20]. They will therefore appear at an higher order in perturbation theory. If we were 
to consider the tree level trispectra {Thhhh, Thhhm, Thhmm, Thmmm), the degeneracies leading to our 
current seven bias parameters get broken and we end up with fourteen bias parameters. It is plausible 
however, given the huge amount of data and highly non-trivial functional form of the trispectra, that 
these fourteen bias parameters can be well determined by just looking at the trispectra (and maybe 
also at the bispectra), so that they can be used for predicting the power spectra and potentially also 
at some level the bispectra. Unfortunately we do not have at our disposal these trispecta data to 
compare with, but it seems however quite clear that, given the huge amount of data, it is possible to 
measure the required bias parameters that we mention from data, without losing much cosmological 
information. The relevant formulas and the estimates for the fc-reach are presented in sec. 2. 

Before passing to the comparison with simulation data, we will extend the description of biased 
tracers in the EFofLSS of [12] to the case where baryons are present and to the case where the 
primordial non-Gaussianities are non-negligible. We will see that interesting new terms are allowed 
in these cases. 

Finally, in sec. 4 we compare our predictions with numerical simulations. We find that the results 
are quite nice, and in agreement with the theoretical expectations. 


2 Biased objects in the EFTofLSS 

2.1 Overdensity of biased objects 

We start from the relation between halos and the underlying dark matter held. As we discussed in 
the introduction, the halo overdensities at a given location depend on the dark matter held and its 
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derivatives, in the combinations allowed by the equivalence principle, through their past trajectory. 
Following [12], we have for the overdensity of halos in Eulerian space: 


4(x,f)~ j dt' H{t') [cs(t,t') :5(xfi,F): (1) 

+Cs2{t,t') : (5(xfi,F)^ : +c«2(f,F) : s^(xfl,F) : 

+C53(t,F) : (5(xfi,F)^ : +C 5 ^ 2 (f,F) : (5(xfi, t')-s^(xfi, t') : +c^{t,t') : '0(xfi,t') : 

+C5st{t,t') : 5(xfi,f')sf(xfi,f') : +Cs3(f,f') :s^(xfl,F): 

+c,(f,F) e(xfl,F) 

+Q 5 (f,t') : e(xfi,F)5(xfi,F) :+c,s(f,f') : es(xfi,f') :+Cet(t,t') :et(xfl,F): 

+c, 25 (t,F) : e(xfi,f')^(5(xfl,F) :+Ce 52 (f,f') : e(xfi,f')(5(xfl,:+Q s 2(^,^') : e(xfl,F)s^(xfl,f') : 

: es(xfi,F)^(xfl,t') : : ef(xfl, F)5(xfl, F) : 

1 

+C 925 (f,F) -^5(xfl,F) + ... , 

Km 

Let us explain the notation. Xfj is dehned recursively as 

Xfi(x, r, r) = X - y dr” v(r", Xfi(x, r, r")) . (2) 

To linear order, one can see that div'^ = —%-d, where "H = aH, and we use the notation' = d/dr, r 
being the conformal time. As usual in Standard Perturbation Theory (SPT), the velocity held can 
be redehned as 6* = such that 9 = 6 ai linear level. For the sake of simplicity, one 

may also redehne (j) so that d'^cj) = 6, we do so henceforth. One can use the equations of motion to 
simplify the number of independent variables, as explained in [20]. Instead of employing n® as the 
independent variable of choice, one can replace 6 with 

h(xfy) = 6'(x,t) - (5(x,t) , (3) 


and djv'^ with Uj, which is dehned as 


fy(x,t) = diVj{^,t) - -6ij9{^,t) - Sij(x,f) . 


( 4 ) 


The quantities rj and will start out as already second order in perturbation theory. We note 
here that, because of the fact that vorticity is generated only at very high order in perturbation 
theory [4, 8], tij can indeed be considered symmetric in i and j up to a high order in perturbations. 
The traceless tidal tensor Sij is dehned as 




( 5 ) 
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It’s easy to show that, to second order in perturbations, 77 ( 2 ) = consequently one can 

introduce a variable 'ijj that is non-zero only starting at cubic order. The expression for reads 

t) = ri{x, t) - ^s2(x, t) + ^^(x, tf . (6) 

Since calculations presented here go up to one-loop order, the quantities introduced this far will 
suffice for all that follows. For later convenience, we introduce non-vanishing products of the following 
variables: 


S^(Xfl,f) = Sy(Xfl,f')s*^(Xfl,f') , S^(Xfl,f) = Sij(Xfi,f')s*'(Xfl,F)sd(Xfl,f') , (7) 

sf(xfi,f) = Sij(xfi,F)F^(xfi,F) , es(xfi,t) = (xfl, F)s*^(xfi, F), ef(xfi,f) = ey(xfi, F)f^(xfi, F) , 

where indices are lowered and raised with (5*h The trace part is absorbed into e, so that eij can be 
treated as traceless. In the last line of Eq. (1) we have the hrst term coming from the bias derivative 
expansion; it corresponds to powers of {k/ku)"^. The scale ku here is the comoving wavenumber 
enclosing the mass of an object [12]. It follows from this that the derivative expansion terms will be 
of higher importance in the case of more massive halos. Finally, the notation ; O : stands for normal 
ordering: 

: 0 := 0 -{ 0 ) . ( 8 ) 

The time-dependence of each order in perturbation theory depends on the order (n) itself and on the 
nature of the held (velocity, overdensity, etc.), but importantly does not depend on the wavenumber k. 
This is due to the fact that the speed of sound for dark matter is only a perturbatively small 
effect, used as a counterterm in the perturbation-theory counting of the EFTofLSS. We can then 
symbolically perform the time integrals in equation (1) to obtain [12]: 


5fe(k,f) = 

= Cyi(t) +Cy2(t) + Cy3(f)5(^)(k,f) + €5,3^^ (k, f) 

+ [C5,l{t) - Cs^2{t)] [(9*5^^^ + [C5,2(^) “ C5,3(f)] 

^Cyi(t) - Cs,2{t) + ^Cs,3{t) 


X 




+ 


+Cs2^i{t) [(5^]k^(t) + Cs2^2{t) [5^]k^(t) - 2 [cs2^i{t) - Cs2^2{t)] 
+c« 2 ,i(t) [s^]^^\t) + c, 2 , 2 (t) [s^]k^(t) - 2 [c, 2 _i(t) - c, 2 ^ 2 {t)] 

+Cst,i{t) [st]k^(t) + Cp,i(t) + C53[(5^][f^(t) -t- C5s2Mk^(t) 


(9) 
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+Cs3 [s^]^ ^ (t) + [e]k + cs [5e] jj. \t) + .... 


We can reorganize this expression pntting together the operators mnltiplying each bias coefficient: 


(Cg(k,t) + Cg(k,t) + Cg(k,t) 
+ cs,2{t)[ Cg(k,t) + Cg(k,t) 

+ Cs^sit) 

+ {t) eg (k, t) 


+ C52,i(t) (cg,(k,t) + cg,(k,t)j + C52,2(t) Cg2(k,t) 

+ C53_i(t) Cg^(k,t) 

+ Cs2 i(t) ^C^zgk, t) H-C^a^k, +Cs2 2(t) cg2(k, t) 

+ c^3,i(t) c3gk,t) 

+ c,t,i(t) C2,\(k,t) 

+ c^,i(t) Cg^(k,t) 

+ C5^2^i(t) Cg^^(k,t) 

+ ca25^i(t) C2g(k,t) 

+ C,(k,t) 

+ C5e,l(t) cg;^(k, t) + . . . , 


,( 3 ) 


,( 3 ) 


( 10 ) 


where the explicit expressions of the Cj operators in terms of the ones in eqnation (9) are given 
in Appendix A, specihcally in eqnation (68). The index n in means the order in pertnrbation 
theory of the operator. Operators describing stochasticity effects are labeled by “e”, we will retnrn 
to disenss them more in detail later on. As pointed ont in [12] not all of these operators contribnte 
in an non-degenerate fnnctional form, and it is nsefnl to check for possible degeneracies 

Note that degeneracies may ocenr only amongst terms of the same pertnrbative order. In order to 
get rid of these rednndancies, we need to hnd the optimal operators basis that spans the fnll fnnctional 
space we have obtained in eqnation (10). One can see that once a prodnet operator is introdneed at a 
given order in pertnrbation theory, it will appear with a different bias coefficient as it is evalnated at 

^It is interesting to observe that one could in principle neglect the study of degeneracies among the operators 
without significantly affecting the performance of the theory against data. This is so because if we use a degenerate 
basis, we span the same functional space as a non-degenerate one. The only aspect that changes between using a 
linearly-dependent basis and an independent one is the actual value of the reduced the number of independent 

degrees of freedom might be overcounted if one uses a degenerate basis. But, once the available number of data is 
much larger than the number of parameters, as is our case, when we compare with several 2- and 3-point functions, 
this is hardly a noticeable effect. 


7 



an higher order in perturbation theory (due to the non-linear couplings or by being multiplied by the 
fluctuating flelds contained in the xg terms). However, as we move to higher orders in perturbation 
theory, new product operators are introduced for the first time. Since the new product operators 
might be degenerate with the lower ones evaluated at higher order, there will necessarily be some 
freedom in choosing the independent operators making up the basis. One convenient choice is to 
keep the bias operators that appear for the first time in the expansion at a given perturbative order 
only if they are not degenerate with any of the bias operators introduced earlier. In this scheme, 
once a specific independent bias operator is introduced, one keeps all of the subsequent ones from 
the same family that appear at higher perturbative orders, unless of course they are degenerate with 
already introduced operators (in other words, we keep the independent ‘descendants’, e.g. is a 

descendant of and it is kept in the basis unless it is degenerate with similar descendants; 

is a descendant of and is kept in the basis unless is degenerate with similar descendants, and 
so on. We dub this choice of operator basis as the ‘Basis of Descendants’ {BoD) one, as, once an 
element is accepted in the basis, its descendants are by default automatically accepted 

Interestingly, it turns out that up to the third perturbative order in the flelds, once the degenera¬ 
cies are dealt with, local in time operators make up the full basis (see also [15]). This is an artifact of 
the lower perturbative orders and will no longer be true at higher orders. In what follows we adopt 
the BoD choice for constructing the operators basis, obtaining: 


(l)st order: 

( 2 ) nd order: 

(3) rd order: 
Stochastic: 


(p(i) 

cfl c 
c®, c 


,(2) 

(p(2) 1 

1 

(5,2 5 

^2 IJ 

r 

,(3) 
(5,2 5 

(p(3) 
^5,3 5 

(p(3) (p(3) (p(3) (p(3) /p(3) 1 

^< 52 , 1 ’ H2,2’ ^3,1’^5,3,, ^*2^2/ 


( 11 ) 


In appendix A, eq. (72), we show explicit linear combinations of the rest of the operators in terms 
of the BoD basis. In practice, it is quite straightforward to check the degeneracy of the operators. 
For example if we are checking the dependence of three second order operators 6 i(ki, k 2 ), 62 (ki, k 2 ), 
03 (ki,k2) we build a linear combination a Oi(ki, k2)-|-/3 02 (ki, k 2 )-l- 63 (ki, k 2 ) and require it to vanish 
for two independent doublets of vectors (ki, k 2 ). This determines a and jS. We then simply ask if the 
linear combination, using the solutions for a and f3 that we just found, vanishes for generic (ki, k 2 ). 
This procedure can then be generalized to combinations of higher order operators. It also applies to 
degeneracies between bias operators and the operators that appear in the effective stress tensor of 
dark matter. For example, there is a degeneracy between the ^ operator, originating from the 

last term in equation (1) and the operator coming form the dark matter held counterterm. Up 

■^An alternative, equally valid, choice could have been to keep the independent operators that appear in the local 
in time approximation (where one ignores the time evolution of bias coefficient) and then supplement them, if needed, 
with a properly chose minimal set of non-local in time operators necessary to complete the basis. 



to cubic level, we are therefore led to write 


( 12 ) 


4(k, t) = (Cg(k, t) + Cg(k, t) + Cg(k, t) 

+ (i) (k, i) + cs,2(3) (i) (k, i) 

+ C5,3 (t) (k, t) + (t) (k, t) 

+ Cs2j(2)(t) + Cs2j(3)(t) C^f2\(k,t) 

+ Cs2^2(t) cja^k, t) + c^ 2^2(^) cj^k, t) 

+ C 53 _i(t) C§yk,t) 

+ Ce4(t) C[^i(k, t) + ... 

where the new bias coefficients Cj are given in terms of the old ones from equation (10) in (73) of 
Appendix A. 

In order to simplify the organization of the computations, it is usually useful to write the pertur¬ 
bative expansion of the dark matter overdensity held in Eulerian framework: 


S(k,t) = 5(^)(k,f) + 5(2)(k,t) + 5(3)(k,t) + 5g)(k,f) + 


where we have 
d(^\k,t) = 
d(^\k,t) = 
S(^Hk,t) = 


d^qi 

(27r)3 




(13) 


(14) 


^fnqnq2)<^2^(k-qi, 

^s^^^^i^^2,^3)dD\k-q3-q2-q3)S^^\q^,t)S^^\q^,t)S^^\q3,t) , 


and Fs^'^ are standard Eulerian perturbation theory kernels (see e.g. [23] for dehnitions). By analogy, 
for the overdensity held of biased tracers it is useful to rewrite (10) as 


(k, t ) = (k, t) + 5;,^^ (k, t ) + 5),^^ (k, f ) + (k, t) + ... + (e-terms) 

where we have introduced the new kernels for biased tracers 

d^\k,t) = j Ki^'>{q^)6^j^\k-q^)6^^'>{q^,t) , 
dh\k,t) = j Kf\q^,q^)5^o\k-qi-q2)d^^\qi^)d^^\^2:t) . 

d qi d q2 d Ki^)(q^^ci2,q3) S^D\k-q^ - q^ - q^) 6^^\q3,t)6^^\q2,t)6^^\q3,t) . 


d2), 


d3)/ 


i:(3) 


(15) 


(16) 


5 f>(k,() = 


(27r)3 (27r)3 (27r)3 
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( 17 ) 


The newly introduced kernels can be expressed in the BoD basis as 

= cyi(t) , 

= cs,i{t) c^ 5 ^iHqi,q 2 ) + £5,2(2) W (I2) + £52,1(2) W 42li(qiiq2) > 

Hqn q2, qs. t) = cs,i{t) q2, qs) + £ 5 , 2 ( 3 )W ^55(qi,q2, qs) + cs,3(t) 42(qn q2, qs) 

+£52,i( 3)(^) 4?i(qi> q2) q3) + £ 52 , 2 ^ c^52l2(qi) q2, q3) 

+£53,i(t) c^s^qi,q2, q3) + 5^2,2(t) c^J^qi, q2, q3) • 

We see that the kernels are linear combinations of bias parameters c„,(i), that are in turn linear 
combinations of the bias c..(f) in eq. (10), as well as subkernels £("'^(5'i, (? 2 , • • •); that are dehned in 
Appendix A (see eq. (70) and eq. (73)). 

2.2 Construction of observables: one-loop power spectrum and tree- 
level bispectrum 

We are interested in computing the set of observables which can later be compared to Wbody 
simulations data. In doing so, we will also extract the values for all the bias parameters. We 
compute the halo-matter cross power spectrum Phm, the halo-halo auto power spectrum Phh at one- 
loop order as well as the matter-matter-halo bispectrum B^mh, matter-halo-halo bispectrum B^hh 
and halo-halo-halo bispectrum B^hh at tree level. 

For the halo-matter cross power spectrum Phm at one-loop order we have: 

PUk,t) = (i“>(k.*)4"(k.())'+ ('5'"’(k.()4"’(k.())'+ ('5“’(k,()4“’(k.«))' 

+ (i <=“) (k, i)4‘’ (k. ()>' + ('5“' (k, ()C. (k. *))' + {€’(k. i)^ ’ (k. i)>' 

= C5,l(t) Puik-,t,t) 

+ 2 / ^^^F®(k-q, q)A:®(k-q, q) Pii(g; f, f) Pn(|k - q|; f, f) 

+ 3 Pii{k;t,t) j (c 5 ,l{t)F^^\]^,-q,q) + Ki^\]^,-q,q)^Pll{q]tp) 

+ (£5,3.,(t) + £5,i(t)) (-(27r)c2(^)(f)) ^Pn{k;t,t) , (18) 

%L 

where the primed brackets {)' means that we are dropped a momentum conserving Dirac 5-function 
from the expectation value. Notice that since the kernels contain unknown bias parameters, 
therefore they need to be expanded. For example, for the convolution term in the second line in the 
last equality in eq. (18), we have: 

2 J PfHk-q,q)A:f)(k-q,q)Pii(g;f,f)Pii(|k-q|;f,f) = 
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Cs,i{t) 2 


d^q 

d^q 


- q, q)c 54 (k - q, q) Ai(g; t, t)Pii(|k - q|; t, t) 


+C52^1(2)(t) 2 


+C5,2(2)(t) 2 j F^‘^\k - q, q)'^^lik - q, q)Pn{q;t,t)Pu{\k ~ q\;t,t) 

77 .( 2 )/,. „ ^H2) 


(27r) 


- q, q)pjyk - q, q)Pii(g; t, t)Pii(|k - q|; t, t), (19) 


A similar expression for the second term is given in eq. (74) and eq (75) in Appendix A. Similarly, 
for the halo-halo auto power spectrum Phh we have: 

Phk{k,t) = (4‘'(k,()4''(k,«))'+ (if (k,()if (k,())' + 2{4'>(k,()4“>(k,())' 

+ 2{ii‘>(k,()ig.(k,())' + 4(|£]q£|,>' 

f 


(2jr)7 


[A:f)(k- q,q)] Pii(g; t,f)Pii(|k - q|; t, t) 


+ 6 cyi Pii(k;f,f) j Kf\k,-q,G^Pii{k]t,t) 

J.2 

+ 2 cyiC 5 , 3 ^^(f) (-(27r)c^(i)(f)) p-Pii(/c;f,t) + Const,, p 

^NL 

d^q 


+ Const,,p J j^Pn{q;t,t) 


( 20 ) 


Notice that the term containing the Ks kernel contributes similarly to both P^^n and Phh and no 

( 2 ) 

new terms need be computed. This is not the case for the Ks ^ kernel where the halo auto power 
spectrum Phh brings new contributions relative to the ones in the cross power spectrum Phm- An 
expression for the second term of the hnal expression in (20) is given in eq. (75) in Appendix A. 
Before proceeding further, let us consider the contributions to the power spectrum arising from auto 
correlations of the stochastic held e as well as potential cross correlations of matter helds and the 
stochastic held e. In case of the halo-matter cross power spectrum Phm, as shown in [12], the leading 
stochastic contribution has a scale dependence: 

, 1^2 

dt' H{t') C,(t,T) (6(k,T) W>EFT-one-loop = (fc/fc3Jl/2 ^ ’ (21) 



where 7 is expected to be an order one number, and kg is the scale connected to the constant stochastic 
contribution which we discuss further below. We can see that this term is highly suppressed relative 
to the rest of the terms at one-loop order in the regime of validity of one-loop results, so we will 
neglect it here. Doing so may no longer be correct at higher loops order, where this and higher 
stochastic terms should be appropriately taken into account. The same term also appears in the case 
of auto power spectrum Phh because of the correlation of the stochastic bias with the stochastic dark 
matter stress tensor, but it is again negligible due to similar arguments (see [ 12 ] for details). 
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Last, we consider the auto correlations of the halo stochastic terms which contain contributions 
of the form Constg ([4 [ 4 )' ~ {{ 27 i)/ks)^. This dehnes a new scale kg which corresponds to the 
mean distance between the biased objects of a given type (e.g. halos of specihc mass). In contrast 
to the aforementioned terms (which scale as k"^), this one generates a constant contribution at all 
scales. We may obtain an estimate for {{27i)/ks)^ reasoning along the following lines: a constant term 
well described by a Poisson-type shot noise contribution will scale like 1/h where n is the number 
density of halos in the specihc mass bin. As mentioned above, 27r/fcs then represent the mean 
distance among halos of the same kind. There exists higher derivative stochastic biases, such as for 


example ( 
is expectec 




[ 4 ) 


(27r)^ 

r .3 

^/T ^ « 


Notice that the scale suppressing the higher derivative corrections 


not to be the inverse distance of the objects kg, but the inverse size of the objects /cm, or 
the non-linear scale /cnl- Plugging in the numbers for our simulation (see sec. 4) we arrive at a value 
of the order ((27r)//cs)^ ~ 10^ — 10^(Mpc/h)^ depending on the bin. Correspondingly, the value of 
the contribution in eq. (21) will be of the order (10^ — 10^)(Mpc/h)^ X {k/k^i^Y which shows why 
we have treated this contribution as negligible. We have verihed that these estimates are consistent 
with our hndings in sec. 4, table 1. 

The bispectrum is given by correlating three helds, each one being either the matter or the halos 
overdensity. For the three combinations of interest we get: 


5mm/i(kl,k2,k3,t) = 2 C5^i{t)\^F^^\]^2,'^‘i)Pl\{k2]t,t)Pii{k:i]t,t) + {k2 -)■ ki)^ 

+ 2 iF^^^(ki, k2, t)Pii{ki] f, t)Pii{k2] t^t) + ... 


— 2 Cyi(f) ^4^^(k2, ^^^)Pll{k2] t, t)Pii(fc3; t, t) + (^2 ki) 

+ 4p,s(ki, k2) Ai(/i;i; C t)Pn{k 2 -, t, t)) 
+ 2cy2(2)(t) '(^ll^{ki,k2)Pii{ki]t,t)Pu{k2]t,t) 

+ 2cs2^^2){t) 4 ?i,s(ki,k 2 )Pii(fci;t,f)Pii(fc 2 ;t,t) , 


^21 1 ^ 3 , t) 


2 c|3(t)P4(k2, k3)Pii(fc2;C t)Pii{k3-, t, t) 


+ cyi(t) ^2 iF4(ki, k2)Pii(/ci; t, t)Pii(fc2; C t) + (^2 ^ 3 )^ + ... 

= 2 cl^{t)(^F^'^\k2,k3)Pn{k2-,t,t)Pn{k3-,t,t) 

-f 4p,s(ki,k2)Pii(fci;t,t)Pii(/c2;t,t) + {k2 fcs)) 

+ 2 Cyi(f)cy 2 ( 2 )(^) (^s,2,s0^^i^‘i)^^^iki,t,t)Pii{k2',t,t) + {k2 —)■ ^ 3 )^ 

+ 2 Cyi(f)c52 ]^('2)(f) ^c^2^]^^5(ki, k2)Pii(fci; t, t)Pii{k2', t, t) + (k2 ^ 3 )^ ) 

4-2c,yi(t)Const,^P Pii(/ci;t,t) , 
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S,j/j,,(ki,k 2 ,k 3 ,t) = 2 K^‘^\'ki,'k 2 )Pii{ki,t,t)Pii{k 2 -,t,t) + 2 permutations 

= 2 c|i(t)^ 4 p,s(ki,k 2 )Pii(/ci;t,t)Pii(fc 2 ;t,t) + 2 permutations^ 

+ 2c|i(t)c5,2(2)(t) (^^ 55 ,^(ki,k 2 )Pii(fci;t,t)Pii(fc 2 ;t,t) +2 permutations) 

+ 2 c|i(t)c 52 ^i( 2 )(t) ,,(ki,k 2 )Pii(fci;t,t)Pii(/c 2 ;t,t) + 2 permutations) 

+ 2cs,i{t)ceS,i{t)Conste,p(^Pii{ki;t,t) + Pii(A:2;t,t) + Pii(/c3; ^, ^)) 

+ Conste^s . (22) 

The bispectrum of the halo-halo-halo density has a stochasticity contribution Constg^s. This quantity 
can be related to Const^^p. As usual with quantities that are renormalziaed in perturbation theory, 
this matching can be done only after the UV behavior is taken into account and renormalization 
performed (for an early discussion of these subtleties in this context, see [2]). We shall do this in the 
following subsection. 

2.3 UV dependence and renormalization of bias coefficients 

In order to correctly take into account the effect that short distance physics has at long distances, 
the bias parameters need to be renormalized. This is arrived at by readjusting the bias parameter 
so as to cancel the contribution from perturbation theory that is obtained from a regime where it 
cannot be trusted, and readjusting it so that the effect of short distance physics at long distance is 
correctly reproduced [20]. The results for dark matter power spectra [3, 9] and bispectra [10, 11] 
seem to show that there is perturbative control up to the relatively UV scale of /c ~ 0.6hMpc“^. We 
will therefore take the point of view that modes with wavenumber higher than that are not under 
perturbative control. At about k ~ 0.6hMpc“^, the matter power spectrum is well described by a 
scaling behavior with slope n ~ —2.1. By approximating our universe with a scaling one with the 
same slope, we can estimate which ones, among the various contributions, are UV sensitive [3]. This 
can be done by Taylor expanding m k the integrands in the former expressions, and identifying 
the ones that are UV divergent if the power spectrum has a slope n ~ —2.1. 

Let us hnd these contributions starting from the terms in eq. (74). We can ignore the loop 
terms that appear already in the calculation of the dark matter correlation functions, as they are 
renormalized by the counterterms present in the equations for dark matter. We therefore can focus 
on the other terms. The only divergent terms are the leading ones in the Taylor expansion and lead 
to the following expression: 


3Pii(A;;f,f) [ A:p)(k,-q, q)Pii(g;f,f) 

J / UV 


2 ,,/ 13 „ 34 „ 47 „ 26 „ 136 „ \ 

^ ^ I ~ 20^2 -|- ^C52 2 -f -^Cs2 2 -f 3 c53 1 1 Pu{k; t, t) 


(23) 
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where 


= ((5(x,t)2) = J Pniq;t,t). (24) 

To appropriately parametrize the UV-sensitivity this expression, we need to readjust one of our bias 
parameters to change this contribution and make it correct. This is possible precisely because the 
form of the UV divergent terms here is the one of the linear bias. We may then renormalize the 
linear bias coefficient so as to absorb it. In practice, one may imagine cyi is a bare parameter, the 
sum of a hnite part and a counterterm: 


^<5,1, finite(^) T counter(^) 5 


(25) 


where counter(^) IS choseu to cancel the contribution from (23), while finite(^) will give the 
correct contribution of short distance physics at large distances. 

After renormalization we have: 


3 Pii{k;t,t) 


d?q 
(27r)3 
d?q 


-q,q)Pii(g;t,f) 


finite 


3 cs,iPii{k] f, t) j (^ 54,5 (k, -q, q) + i) Pu{q] t, t) 


+3 (cya + 15cs2^2)Pii{k-,t,t) 


M3) n t 47 


(277) 


^ 5 , 3 ,.(k, -q, q) - ^ ) Ai(g; t, t) 


(26) 


Here with the subscript gnite we mean that we have removed the UV dependent part. This ‘hnite’ term 
is however meaningful only when added to the counterterm. In other words, only loop terms plus 
counterterms give sensible results that are independent of the UV physics and the renormalization 
procedure. From (26) we see that renormalization brings to light one additional degeneracy between 
the bias parameters £ 5^3 and £^2 2 , which appears both in cross and and auto power spectrum. Similar 
behavior appears in convolution-type contributions of the halo-halo power spectrum in equation ( 20 ). 
As an example, we can look at the term proportional to the £^ 21 ( 2 ) coefficient, which is of the form 


where 




(Pq 




UV 



(27) 


(28) 


In the UV these contributions will be analytic functions of k and as such they can be reabsorbed 
in the dehnition of the stochastic coefficient Constg^p and the higher derivative ones (at the order 
at which we are working, there is no need of the higher derivative ones). Similarly to eq. (25), the 
latter will then be 


Constg^p Constg^p gnite T Constg^p counter 7 


(29) 
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where Conste^p counter is chosen to cancel the contribution from eq. (27). Note that, by symmetry 
reasons, Phm does not have the stochastic Constg^p contribution. Indeed all the convolution-type 
divergencies cancel exactly so the hnal result is hnite. 

Note that only the hnite part of Constg^p corresponds to the physical contribution and is thus 
connected to the stochastic contribution in the bispectrum, Conste^p. If, as it is expected, we assume 
a Poissonian nature for the statistics corresponding to the mean halo number density, we have 
Constc,p = Const^p , which is the well known relation of the shot noise contribution in the 
power spectrum and bispectrum (see e.g. [21, 40, 22]). Equipped with this relation, one is able to 
reduce the number of stochastic bias parameters down to one. 

At this stage we can give the hnal set of bias parameters that describe all the introduced statistics. 
We perform the following replacements: 


cs,i -t bs,i - cr^(t) 

Cs,2 (2) —^ ^( 5,2 , 
cys + 155^2^2 bs,3 , 


34 


Cs,2 + 


47. 26. 

+ 2C521 -f —C52 2 + 


136 


Co2 


s^,2 


-p 3 c^3 ] 




Cs,3cs 5 

CeS,l b^S,l , 

Conste^P —>■ Conste^l — 

-E{tr (~ci+ci 2(2) + ^%,1(2) ~ 2 C5^lCy2(2) ~ 2 CyiC 52 ^i( 2 ) -|- 2 Cy2 (2)C52^i ( 2 ) 


(30) 


where the bias coefficients bs^i, 65 , 2 , ^ 6 , 3 , bs^, 652 and Constg are hnite contributions and we have 
used the relation Const^^p = Const^ The cross power spectra of biased tracers and dark matter is 
therefore given by: 


Phm{k,t) =65,i(t) (^Pu{k;t,t) + 2 j Ff ^(k - q, q) - q, q) Pii(g; t, t)Pii(|k - q|; f, f) 

(Pq 


+3 Pii{k] f, t) 


( 2 vr) 


,-q,q) + c} /,(k,-q,q) + Pii(g;f,f) 


(Pq 


hAt) 2 J Pf)(k-q,q) (^Pf)(k-q,q) -44,,(k-q,q)) Pii(g; f, f)Pii(|k - q|; f, f) 

+ bsAt) 3 Pii(fc; k t) J “ i) 


^Notice that strictly speaking, in the renormalization of bs^i and Const one could introduce additional contributions, 
such as those arising from operators such as Se^-^, s^e, 6^e, .... However, all of these terms do not contribute to the 
bispectrum at the order we are working on, and are trivially reabsorbable by a shift in Const and bs^i- We therefore 
omit to write them explicitly. 
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( 31 ) 


+ bs2{t)2 j ^^^Ff)(k-q,q)Pii(g;t,t)Pii(|k-q|;t,t) 
+ {h^Xt) - 2(27r)c2(i)(t)65,i(t)) ^Pii{k]tX) • 

%L 

Note that the following relation holds: 



<Pq 

{2nf 



,(k,-q,q) + 


c(3) 

^S,3,s 


(k,-q,q) - §) Pii{q;t,t) 


and for {65^1 —)■ 1, bs,2 —^ bs^i, 65,3 —)■ &5,i, bs2 —)■ 0, bc^ —)■ 0} our formula reduces to that of the 

standard dark matter case (below at ( 62 )), as it should. 

Similarly, for the auto power spectra of biased tracers we have: 


Phhik,t) =bliit) Pii{k;t,t) + 2 


d^q 

(27r)3 


4p,.(k - q, q) ^11 (g; t, t)Pii(|k - q|; t, t) 

d?q 


+6 Pii(A;;t,t) J + l) 

&5,i(t)&5,3(^) 6 Pii(k;t,t) j (^c^55.s(k,-q,q) - i) Ai(g;t,t) 

Hk - q, q) - 4p,s(k - q, q) Pniq;t, t)Pii(|k - q|; t, t) 


k‘s,2{t) 2 

blit) 2 


d^q 

(27r)3 

d?q 

(27r)3 


Pii(g;t,t)Pii(|k- q|;t,f) 


b&,i{t)bs^2{t) 4 


d^q 


( 2 vr: 


p( 2 ) 

,3 ^<5: 


l.(k-q,q) (PfHk-q,q) 


q,q)) Pii(g;t,t)Pii(|k 


q|;t,t) 


+ &5,i(t)&52(t) 4 J ^^^41s(k-q, q)^ii(g;^,f)Ai(|k-q|;t,t) 

+ fe5,2(t)fe52(t) 4y (^pf)(k-q,q)-44,,(k-q,q)) Pii(g;t,t)Pii(|k 
k^ 

+ bs,i{t)bs,3^^ (t) 2 p-Pii(A:; t, t) + Const, 

— S(t)^ + 6^2 + ^52 ~ 2 ~ 2 65,1652 + 2 65,2652) . 


q|;t,t) 


( 32 ) 


Again, as a check, taking bias parameters {65,1 —?• 1 , 65,2 —)■ 65,1, 65,3 —)■ 65,1, 652 —)■ 0 , 6,,^ —)■ 0 } we 

reproduce the dark matter power spectrum result ( 62 ) below. 

Let us consider the effects of renormalisation on the bispectrum expressions in eq. ( 22 ). Here too, 
in analogous way to eq. ( 30 ), we employ renormalized coefficients. Notice however that we have to 
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include only the finite part of the coefficients (schematically shown in eq. (25)). We do so because 
the bispectrum terms are in no need of being renormalized at tree level. The ‘counter’ part of eq. 
(25) is present in order to correct for the incorrect contribution from loop diagrams. Since we do not 
include loop diagrams in the bispectrum, we should not include the ‘counter’ part of the counterterm, 
as at this order, there is nothing to correct for. This is the usual fact that tree-level expressions do 
not need to be renormalized. In summary, and schematically, at tree level for the bispectrum we 
perform the following replacement (for cyi and correspondingly for the rest): Cyi —)■ finite = 


2.4 IR-resummation 

As detailed in [7, 12], in the EFTofLSS, after the renormalization is performed, the loop expansion, 
complemented by the insertion of the relevant higher order bias coefficients, corresponds to an ex¬ 
pansion in the parameters that regulate the dark matter expansion, those are e 5 < and es> [7, 9]. 
These quantities are defined according to 


6 s> h 


^5< = 


d^k' Pii(fc') 

(27r)^ 


d^k' 

( 2 vr)= 




(33) 


where Pii {k) stands for the dark matter power spectrum. For a given wavenumber A:, es> will 
represent the effect of the displacement due to short wavelength modes, while es< describes the effect 
of tidal forces due to long wavelength modes. Both these quantities scale proportionally to /c/Zcnl- 
The bias derivative expansion corresponds to one in powers of {k/kuY, where ku is the comoving 
wavenumber enclosing the mass of an object [12]. Within an Fulerian treatment we ought to expand 
also in the displacement due to long wavelength modes es< = (fc(5s<)^, defined as 


d^k' Pii(fc') 
0 (27r)^ 


(34) 


As explained in [9], in the fc-range of interest es< is of order one; as a consequence one cannot Taylor 
expand in this parameter. Fortunately, for equal-time, dark-matter correlators the only effect of es<_ 
is to affect the Baryon Acoustic Oscillation (BAO) peak. In Fourier space, this shows up in about 
2% residual oscillations that are not correctly reproduced if one Taylor expands in es<. In this work, 
we aim at exploring the UV reach of the FFTofLSS for biased tracers, and therefore we perform the 
IR-resummation that corresponds to treating es< non-perturbatively only for the power spectrum. 
In this case, since there is no velocity bias that is not higher derivative [12, 15], the formulas for 
IR-resummation simply extend the ones for dark matter [9, 12]. 


P&aSbik] ^1,^2) I TV 



M\\n-j (k, k'] ti, fa) Ps^Sbjik'; fi, ^2) 


(35) 
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( 36 ) 


where M\\^_.{k,k']ti,t2) is defined as the double Fourier transform of -Pint||iv-i(^l'?j^2) 

= j d^q (q. -k; h, fa) , 

and PsaShj stands for the j-th order term for the power spectrum in the Eulerian calculation. We 
refer the interested reader to [ 12 ] for further details on the notation. 

2.5 Estimated theory errors 

An important property of the EFTofLSS, and of every theory that is not a model, is that it allows 
for a reliable estimate of the error at every given perturbative order. We can readily do this by 
assuming a simplified scaling universe behavior corresponding to a slope n ~ —1.7 for the linear 
power spectrum. There are several contributions that come at the next order both in the power 
spectra and in the bispectra, which scale differently accordingly to size of the bias parameters and of 
the stochastic terms. This means that, when this is relevant, for some quantities we quote estimates 
for the contribution from more than one diagram, and we take the largest among those as our 
theoretical error. For the matter power spectrum we thus obtain 

/ 7 \ 2.6 

~ (27r)2 — Pn(fc) , (37) 

V^NL/ 

which comes two-loops of dark matter self-interactions. For the halo-matter power spectrum, we 
take 

APftm ~ (27r)^ Max [ 65 , 1 , 652 ^ 1 ] ( - — ) Pu{k) , (38) 

V^L/ 

which again come from a two-loop diagram using dark matter self-interactions. The first term comes 
from using the linear bias, while the second comes from the quadratic one. For the halo-halo power 
spectrum, we have similarly 

AP,^;, ~ (27r)2 Max [ 6 |i, 6 ^ 2 , 1 ] Pn{k) , (39) 


depending if we use the linear or the quadratic bias. 

Similarly for the tree-level halo-halo-halo, and halo-matter-matter bispectra, we derive the fol¬ 
lowing error estimates, associated to dark matter loops: 


APmmm ~ (2^) 



1.3 


B 


tree level ? 


(40) 
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^Bhmm ~ (27r) Max [ 65 , 1 , 652 ^ 1 ] 


NL 


1.3 


B 


tree level 


(41) 


where k is the largest wavenumber of a given triangle. For tree-level halo-halo-matter, the contribu¬ 
tion from the stochastic bias in eif can become relevant. We therefore write 


AR 


hhm 


Max 


(2?r)Max ,] (2ir) 6,ji (^) 


.(42) 


where fcdm is the wavenumber associated to the dark matter held. Finally, for the tree-level halo- 
halo-halo bispectrum, we take 


AR 


hhh 


Max 


Max[(27r)65_i^, 


1.3 


27r\ V k '' 


, ) .Streelevel) &(5,l^e(5,l ( , II; 

^NL / \^S J \ ^NL 


Pii{k) 


(43) 


Strictly speaking, the relation given just above for the bispectum error is obtained for the equilateral 
triangles conhguration where all three momentum vectors have the same amplitude. Arbitrary tri¬ 
angle conhgurations would make the relation more cumbersome, but without signihcantly changing 
the result. Similarly, for the porpuse of estimates, we work in the local in time approximation, and 
therefore we do not distinguish between bs,i, ^ 5 , 2 , &< 5 , 3 , • • • and we do not take into account of possible 
degeneracies. Importantly, we notice that tracers with high values of their bias coefficients will tend 
to be less well behaved in perturbation theory, in a way that can still be predicted in the EFTofLSS. 
These estimates will be later used to estimate the fc-reach of our calculations with the EFTofLSS. 


3 Baryons and Primordial non-Gaussianities 

The results presented in the former section would allow us to proceed to compare the EFTofLSS with 
data of biased tracers. We delay this to hrst extend the formalism to account for two interesting 
effects, even though we will not study them in the comparison with data in this paper. The derivation 
of bias results in the EFTofLSS, as described in sec. 2, has been relying on a few essential assumptions: 
the validity of general relativity, the presence of only one species in the universe (a dark matter 
component), and the absence of primordial non-Gaussianities. The fact that general relativity holds 
was assumed when, for example, we did not allow the halos in Eulerian space to depend directly on 
the gravitational potential <F or its hrst derivative c?<F. The presence of only one species was assumed 
in the fact that we wrote only one matter component, and that the halos were not allowed to depend 
on the velocity of the particles, apart for the streaming terms. Finally, the assumption of Gaussian 
initial conditions was implied in the fact that the only effect of long modes on the formation of halos 
was simply due to evolution effects, such as the fact that long tidal forces or velocity gradients might 
enhance or delay the collapse of a certain short region into halos. The initial conditions of short 
modes that end up collapsing were not affected by the long modes; this is where the hypothesis of 
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gaussian initial conditions, which indeed implies the fact that modes of different wavelengths do not 
affect each other, has entered. 

In this section, we will generalize the effective description of biased tracers presented in sec. 2 to 
the case in which we include the presence of baryons and of primordial non-Gaussianities. We will 
leave the inclusion of deviations from general relativity to a later time. We will simply provide the 
relevant equations, and we will compare to data in subsequent studies. 


3.1 Including Baryons 

The inclusion of baryon effects in the EFTofLSS was recently performed in [14]. There, it was shown 
that baryons can be described at large distances, even in the presence of star formation events, as 
a second species on top of the dark matter, where the effect of short distance non-linearities and 
of star formation physics at long distances is encoded in an effective stress tensor and an effective 
force. The effective force, absent in the case of only one species, takes into account the momentum 
exchange at short distances between dark matter and baryons through gravity. By comparison with 
measurements of the power spectrum from simulations, it was shown that the EFTofLSS can repro¬ 
duce baryon effects at percent level with a one-loop calculation up to the relatively high wavenumber 
k ~ 0.6hMpc“^, consistent with the estimates of the theoretical errors. We now proceed to the 
theory for biased tracers. 

The generalization of the formulas of sec. 2 is rather straightforward. We start from Eulerian 
space, where the derivation follows in analogy the writing of the effective stress tensor and force 
terms in [14]. In the EFTofLSS with two species, the degrees of freedom that we are able to use are 
the gravitational field, the density and velocity gradients of the two species, and, finally, the relative 
velocity of the two species in the center of mass frame. 

The contribution of long wavelength perturbations of dark matter and of baryons is expected 
to be weighted by the relative abundances of baryons and dark matter: Wh = , Wc = 

This is important because, as shown explicitly in [14], the factor of Wb suppresses perturbatively the 
contribution of baryons. When introducing baryons, there is therefore a sort of doubling of the bias 
parameters that represent the effect of matter fluctuations, even though those associated to baryons 
are perturbatively less important. 

As anticipated, in the presence of two species, the bias can depend directly on the relative velocity, 
which is not zero in the center of mass frame [14, 15]. In fact, in the center of mass frame we have 

^^c,CM = Vc- {WcVc + WbVb) = Wb{Vc - Vb) , Ufo.CM = Vb - (Wc^’c + WbVb) = Wc.{Vb - Vc) ■ (44) 


The halo fields should therefore be allowed to depend on them. Due to the units of the velocity, 
the only combination that involves no derivatives on the velocity field directly (as these terms are 
allowed in the case of single species as well) are quadratic and effectively proportional to Wb, such as 
Wa vl / H, where a, a' run over c, b. 

Summarizing, if in the only-dark-matter case we write, as pointed out in earlier sections from [12], 


4(x,f) ~ 


dt' H{t') 




9^(/)(Xfi, t') 

Hit'Y 


(45) 
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+ ... 


'^^divAt^t ) H{t') ^didj<l>d'^di<j>\ti't ) H{t'Y 

+c,(t,t') e(xfl,t') + c^a^^{t,t') 6(xfl, t') ^ ^ + • • • 

in the presence of baryons we have the following generalization 

aV(Xfi, t') , - . , . . X 

+ ^ + C4(t,f ) Wbd)b(Xflfe) 


-l-f\ 


4 (x,t)~ y dt'H{t') Ca 2 ^{t,t') 

(t, t) w, —-+ ca,yi(t, t) Wb —— 

+caiaj(t>d^a2<f>[t,t) H{t'y ’ 

+c,^(t,t') w;cec(xfi,c,t') + c,^(t,t') wj, eb(xfl,b, ^') 


+c, 


^a^ii^t') w;cec(xfl,c,t')^^^^^ + c.,a2<^(t,t') w;, eb(xfl,b, 


I - u +'^ ^^‘/’(xfl,^') , >y;;^ i I ^ 

+cav(t,t)^ i7(t0^ ^ 2^ w^<7'«^.,CM(xfl.>^ 




cr,cr'=6,C 


if 


(46) 


where 

nr nr 

Xfl,b(x, r, r') = X - / dr" Vfe(r", Xfl(x, r, r")) , Xfl,e(x, r, r') = x - / dr" Vc(r", Xfl(x, r, r")) . 
Jr' J t' 

(47) 

and where the xg term in ^ ^ should be meant as by using the equations of motion 

^ ^ '5c(xfl,c, t) + Wb Sb{xfi,h, t) . (48) 

Notice that we have inserted a factor of Wb in front of the baryonic fluctuations. This is because the 
effect of baryons on halos is expected to be suppressed by Wb- The suppression in Wb is probably 
absent in the case the collapsed objects are galaxies, so if one were to use these formulas for galaxies, 
one should allow for these terms to be non parametrically suppressed by Wb- 

At this point, we proceed in the analogous way as we did for the case of dark-matter-only [ 12 ], 
We Taylor expand in all perturbations, and we perform symbolically all the time integrals to obtain 
independent bias coefficients for each order in perturbation theory; one then finally explores the 
degeneracy between the functional terms, along the same lines we elaborated upon in sec. 2 . 

In the case of two species, the expression for the biased tracers written in Eulerian space, as the 
ones above, expand not only in es<, es> and es<, but also in the parameters associated to the relative 
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displacements [14]: 


'OO 


6 


rel 

s< 


{k) = k'^ 



d^k' Pii(fc') 
(27r)^ 


et{k) = k^ 


d?k’ Pii(P) 

(27r)^ 


(49) 


where Pi\{k) is the power spectrum of the log-derivative with respect to the scale factor a of the 
difference in the dark matter and baryon overdensities: d{5c — 5^)/d\oga. The parameters esc ^md 
are order one in the current universe for the fc’s of interest is order one only at high redshifts), 
and so cannot be treated perturbatively. Fortunately, they are IR dominated, and can therefore be 
resummed. Formulas for the IR-resummation in the case of two species were provided in [9] for the 
baryon and dark matter helds. As the usual IR-resummation for single species extends trivially from 
overdensities to tracers [9], so it does for the case of two species [14], For the overdensity power 
spectrum, we have 

TV 

P“ft*)liv = 5Z / , (80) 

j=0 

where a = hh or hm. The quantity P^{P]t) is the Eulerian result at j-th order obtained by 
expanding in esc, ^s>, ^s<, treated as equal, P°‘{k;t)\N is the IR-resummed result up to N- 

th order by expanding only in esc, ^s> and and treating them as equal. The matrix 
is dehned in eq. (4.5) of [14]. It is important to keep in mind that this formula is obtained by 
treating at non-linear level the displacements of the dark matter and baryon particles, and not of the 
adiabatic and isocurvature modes. It is in this basis that the IR-resummation takes a particularly 
simple form. Notice hnally that for this formulas to be consistent, it is crucial that there is no linear, 
non-derivatively suppressed, velocity bias, as it is indeed the case. 


3.2 Including primordial non-Gaussianities 

In the case when the primordial fluctuations are non-Gaussian, the formulas of sec. 2 are incomplete. 
In fact, in the derivation of those formulas, it was assumed that the only way the long modes affect the 
number density of collapsed objects is by affecting the evolution of the short modes that eventually 
collapse. However, in the presence of non-Gaussianities that correlate modes of different wavelengths 
in the initial conditions, the presence of long wavelength fluctuations can in principle affect the actual 
distribution of the short modes before horizon re-entry, when all evolution effects are still null. The 
explanation of this effect is quite intuitive: suppose that the presence of a non-Gaussian correlation is 
such that, in the presence of a positive primordial fluctuation, the variance of the short fluctuations 
is enhanced. Then, it is quite intuitive that it would be particularly easy for collapse to happen in 
that region. This effect is particularly interesting because the initial conditions are, contrary to the 
evolution effects, non-local in space, and can therefore have a different /c-dependence than the one 
obtained from standard biases 

®The fact that we need a way to generate non-local in space initial conditions is what is named as the horizon 
problem in cosmology and leads to the need of inflation. It is quite interesting to see the resurgence of such an historical 
problem in this unusual setting. 
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The description of bias in the EFTofLSS is based on an hierarchy of scales between the long 
modes ki that affects the collapse of short modes with wavenumber fcg of the order of the one whose 
wavelength encloses a region with the same mass M as the collapsed object. It is therefore an 
expansion in k^/ks -C 1. Since there is no hierarchy in time, the EFTofLSS for biased tracers is 
non-local in time [12, 4, 6], so there is no rigorously justihed analogous expansion in lul/ujs- The 
limit ki/ks -C 1 implies that the functional dependence of the number density of collapsed objects 
on the long-wavelength ones will probe the squeezed limit of non-Gaussian correlation functions. 

In summary, two ingredients are necessary to include the effect of primordial non-Gaussian cor¬ 
relations into the bias. The first is that these correlations are present in the initial conditions, and 
not developed by the time evolution when modes are sub-Hubble. Second, the effect can be well 
described by the squeezed limit, k^ of correlation functions. 

Let us follow through these considerations and derive an equation for biased tracers. We will 
focus our discussion on the case where there is a non-trivial three-point function, the so-called /nl 
case. The study can be straightforwardly extended to the case of primordial non-Gaussianities for 
higher Wpoint functions. In the squeezed limit k^ kg, we can write any three-point function as 
the one induced by the held ® 

C(ks) ~ C,(ks) + /nl (Lf C 5 (ks - ki)Q(ki) . (51) 


where C, is the curvature perturbation that is constant on super-Hubble scales '. (g is an auxiliary 
gaussian variable, a is some power that depends on the particular inhationary model. For example, 
in the local case, f^l, a = 0, for the equilateral, [43], and orthogonal, [44], cases, a = 2 

quasi-single held inhation [45, 46], gives 0 < a < 3/2. Strongly coupled conformal sectors can give 
rise to 0 < a < 2 [47]. If precision requires, additional subleading terms in the squeeze limit can be 
added in (51) in a similar fashion. 

Well after horizon re-rentry, but still early enough to neglect all gravitational non-linearities, the 
primordial density huctuation can therefore be written as 


iolks) +/NL/(ki,iia)i5j(ks - ki,i,„) , 


(52) 


where 


•} t[n) 


5g(kL,tin) 


(53) 


2D{U^) klT{k) 

and where T{k) is the transfer function. The super-script T) signals that this is the initial density 
held that will then be non-linearly evolved up to the present time. It is non-Gaussian simply because 
of the primordial non-Gaussianities (see for example [48]). 


®This formula does not account for anisotropic squeezed limits that are generated in some inflationary models, 
which however can be included using the same formalism as developed here (see [50, 51, 52] for subsequent works). 

^For the proof of the constancy of C at quantum level see [42]. 

^Concerning the squeezed limit on the orthogonal case, there is some confusion in the literature due to the fact 
that in ref. [44] there are two templates for the three point function that is induced by The one in the main 

text has an incorrect squeezed limit, but it is used for CMB analysis where the signal in the squeezed limit is limited. 
The one in the appendix has the correct squeezed limit, leading to a = 2, and should be used for the bias. 
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In general, the number of collapsed objects at a given location, n/i(x, t), will be a function of the 
evolution of the short fluctuations evaluated on the path that led them to the hnal point. The short 
fluctuations will be well described by the set of all their correlation functions. Symbolically, 

nh{^,t) = j dt' K{t,t') fh ...) . (54) 


The subscript ()/ means that the average is taken over the short modes by keeping hxed the long 
modes. Obviously, it is impossible to compute perturbatively this expression. But, as for the case 
of no non-Gaussianities, we can still Taylor expand in the long wavelength fluctuations. Due to the 
primordial non-Gaussianities, the dependence on the long modes is not just the one induced by the 
gravitational dynamics. For example, the two-point function of short modes contains terms that 
scale as: 

(^s)z(Xin,tin) D ((5^)o(tin)/NL0(Xin,tin) • (55) 

This term allows us to deduce that the overdensity of collapsed objects will depend on terms of the 
form 


5fe(x,f) D J (iFiF(f,T)&(T)/NL0(Xfl(f,fin)) = /nL&i 0(Xfl(t,hn)) 
- b\ /nl 


~ ~ B^ B^ 

0(x,fin) - ai 0 (x,fi„)— 0(^)(x,f) + 


where in the last line we have expanded in xg to cubic order. 

As said, the correlation functions of the short modes cannot be computed explicitly. However, 
the correlation functions that are present in the initial conditions will be enhanced by those that 
come from gravitational interactions. Goncerning the effect of the long modes, this means that on 
top of depending on 0(xfl(f, fm), hn), the collapsed objects will depend on the held that were allowed 
to be present in the absence of primordial non-Gaussianities, evaluated on their past trajectoreis, 
plus all products of these helds evaluated on their past trajectories with powers and derivatives 
of 0(xfl(f, fin), tin)- The peculiarity of these terms is associated to the presence of the transfer function 
in (53), as well as potential non-analytic powers of k^: this non-analytic behavior is non-degenerate 
with terms that we could have written as arising from local evolution. Indeed, local dynamics is 
mirrored in analyticity in Fourier space, and therefore any non-analytic behavior must come from 
the initial conditions. 

In the presence of primordial non-Gaussianities, we can therefore write the following additional 
terms on top of the ones in (9): 


Shi^iyt) ~ /nl 0(Xfl(f,fi„),fi„) 


dt' H{t') 




a^ 0 (Xfl,f') 

Hit'Y 


(56) 
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+ ... 






H{t') 


+ c 


didjcfid^d^ <t> 


{t,t 


didj(p{xa,t') 


H{f) 


Hit') 


+c/(t, t') e(xfl, t') + c^Q 2 ^it, t') e(xfl, t') 




I . 0 /. ./N 5V(xfl, t') 


Hit') 


+ 


+ /nL^ 0(Xfl(t,tin),tin)^ / dt'H{t') 




c"(t,0 + CaVt,0 


Hit') 


+ ... 


As in the former sections, the time integral can be symbolically done, to lead to the following new 
bias coefficients: 


d' 


dhik.t) = f^L ( cf(t) [(j)]kit) - cf(t) [di(j)-^e^\it) + ... 


(57) 


&' 


+cyi(^) [^6^%it) - Cg^^it) + cs^2i^) [^6^%it) + 


+ 


- Cs^2it) ^^^%it) + • • • + [^s^]^^\t) + ... 


d2 


■ d' 


+cj,i(t) + • • • + c,)\(t)[0e(^)]fc(t) - cj^iit) [di^-^0^^'^ e^%it) 

+/ nl ^ (cfit) [^%it) + ■ ■ ■ + eg it) [02 6^%it) + ... cgm\^%it) + .. 


where the dots stand for terms higher than cubic order and we have expanded the xg dependence. 

In the local in time approximation, the linear bias coefficients due to primordial non-Gaussianities 
have been known for some time. The linear scale-dependent bias induced by was hrst discovered 
in [37], even though the presence of the transfer function factor was introduced only in [38]. It has 
since been studied extensively in the literature. The presence of the transfer function factor that 
makes the functional form of the linear bias distinguishable from the standard one was pointed out, 
for the equilateral and orthogonal cases, in [39, 41]. We point out here the generalization to the 
non-linear biases, to the non-locality in time, to higher derivative terms as well as, next, to the 
IR-resummation. 

In fact, so far this formulas are written in Eulerian space, which expands in es<. This parameter 
is order one in the current universe, and so needs to be resummed. However, the same considerations 
that apply in the absence of primordial non-Gaussianities apply here: only long wavelength displace¬ 
ments need to be resummed, and, as explained in [12], the IR-resummation for biased tracers is the 
same as for dark matter as developed in [9] 

®The reader concerned by the fact that the variable 0 appears only evaluated at tin, rnight notice that this can be 
thought of as a non-local in time bias with support as a i5-function at an early times. Therefore, the same discussion 
as in [12] applies here. 
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4 Results and Comparison to A^-body simulations 

In this section we compare our results for halo cross and auto power spectrum and bispectrum of 
dark matter halos with Gaussian initial conditions, against iV-body simulations. We £t for seven 
independent bias parameters; six enter in the halo-halo auto power spectrum, hve enter in the halo- 
matter cross power spectrum, four enter in the halo-halo-halo auto bispectrum, and three in the 
halo-halo-matter and halo-matter-matter cross bispectra. We fit these parameters by minimizing 
the joined function. We include below the results and corresponding plots obtained using these 
values. 

4.1 A-body simulations 

We compare our theoretical results with the Millennium-XXL (MXXL) iV-body simulation [26]; it 
evolved 6720^ particles inside a comoving cubical region whose side measures = 3000Mpc/h. 
The evolution onset corresponds to z = 63 and lasts until the present day. The MXXL simulation was 
performed with the following settings for cosmological parameters: a matter density of 12^ = 0.25 
in units of the critical density, a cosmological constant with Ga = 0.25, a Hubble constant h = 0.73, 
a spectral index Ug = 1 and a normalization parameter cxs = 0.9 for the primordial linear density 
power spectrum. The initial conditions were set using the linear theory power spectrum obtained 
from the Boltzmann code CAMB [27]. 

In this work we focus on the auto and cross power spectrum and bispectrum of the halo and 
matter density contrast at redshift z = 0. For the dark matter component, as described in [10], 
we first construct a density held by assigning simulation particles onto a 1024^ cubic grid using a 
“Cloud-in-Cells” deposit scheme. This is followed by Fast Fourier transforming the given density 
held and correcting for the ehects of the assignment scheme by dividing each Fourier mode by the 
Fourier transform of a cubical top-hat window function. 

Dark matter halos are identified using the “Friends-of-Friends” algorithm. Linking length is set 
to 0.2 times the mean particle separation and only halos with more than 20 particles are used in 
the analysis. Since we are interested in how clustering of halos depends on the halo mass, the dark 
matter halo sample is divided into four subsamples (binO, 1, 2 and 3) according to their mass. Mass 
bins are dehned as in [25, 24]. Because of the fact that the scaling of measurement of the bispectrum 
goes as the sixth power of the number of grid cells per dimension, cross and auto bispectrum are 
measured up to modes k < 0.15hMpc“^. 

The cross and auto power spectra are computed by spherically averaging the amplitude of Fourier 
modes in a range of radius Ak. All the bispectra B{ki, k2,0), as detailed in [10] for the dark matter 
case, are obtained by performing a nested loop over all grid points (1024®) and averaging the quantity 
S{ki) 6 {k 2 )S{k 3 ) over triangles whose sides satisfy: 

k{ = ki± ^Ak, ki = k2± = e± ^Ae , (58) 
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Figure 1: Dark matter power spectrum at redshift z = 0 divided by the nonlinear simulation data. The 
one-loop EFTofLSS result is shown after (solid blue line) and before (dashed blue line) performing the IR 
resummation procedure. For comparison, we also show the one-loop SPT result (solid orange line). The 
parameter is determined by fitting the IR-resummed EFT predictions over the range 0.15/iMpc“^ < 
k < 0 . 3/1 Mpc“^. The thin gray dashed lines signal the 2% systematics error assumed for the simulation 
data. 


and where = cos ^ ■ ^ 2 ^ • For the purposes of this work we have set 

Ak = 2Ti/Lhox = 0.0021Mpc//i, SO = 7 r /20 . (59) 

Finally, we review the estimation of the errors in the measurement of the power spectra and bispectra. 
For the cross and auto power spectrum we have for the following: 

= (60) 

where [xy) indices may stand for both, halo-matter (hm) or halo-halo {hh) power spectrum, and n 

is the number of modes contributing to the given /c—bin. Similarly, for the auto and cross bispectra 
one has (see e.g. [29, 28, 31, 30]): 

(fcs) , (61) 

^triangles 

where x and y encompass the halo h or matter m case. Note that above S 123 is equal to 6 , 2,1 
for equilateral, isosceles and general triangles respectively and that the relation = k'l + k^ — 
2 kik2 cos( 6 ') is satisfied. The quantity utnangies stands for the number of triangles contributing to a 
given conhguration, which, in the case at hand, is directly counted inside the bispectrum code. 

Even though derived in the Gaussian limit, these estimates give a better approximation for the 
scales we shall be concerned with if we keep P^vl in heu of Pn. In addition to the error estimates 
presented above, we add a 2 % error for each bispectrum data point, in order to account for possible 
systematic errors in the simulation. 
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0: k 

min - 0 

.04h/Mpc, 

kmax - 0 - 

16h/Mpc 
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hh 

hmm 

hhm 

hhh 


P 

+ 
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— 

0.0272 

1.000 
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+ 
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— 
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+ 
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— 

+ 

— 
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0.9987 

+ 
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— 

— 

+ 
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+ 

+ 

+ 

+ 

— 
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+ 

+ 

— 

+ 
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+ 

+ 

— 

+ 

+ 
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0.2166 
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+ 

+ 

+ 

+ 
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— 

0.0372 

1.000 

+ 

+ 

+ 

— 

— 

0.662 

0.9937 

+ 

+ 

— 

+ 

— 

0.615 

0.9982 

+ 

+ 

— 

— 

+ 

0.730 

0.9724 

+ 

+ 

+ 

+ 

— 

0.849 

0.8911 

+ 

+ 

+ 

— 

+ 

0.846 

0.8963 

+ 

+ 

— 

+ 

+ 

1.17 

0.09115 

+ 

+ 

+ 

+ 

+ 

1.13 

0.1105 


bin 2: kmin=0.04h/Mpc, kmax=0.lOh/Mpc 

hm 

hh 
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hhm 

hhh 


P 

+ 

+ 

— 

— 

— 

0.0232 

1.000 

+ 

+ 

+ 

— 

— 

0.412 

1.000 

+ 

+ 

— 

+ 

— 

0.476 

1.000 

+ 

+ 

— 

— 

+ 

0.504 

0.9999 

+ 

+ 

+ 

+ 

— 

0.627 

0.9992 

+ 

+ 

+ 

— 

+ 

0.640 

0.9988 

+ 

+ 

— 

+ 

+ 

1.04 

0.3656 

+ 

+ 

+ 

+ 

+ 

1.02 

0.4101 


Figure 2: p-values are given for best fit bias parameters procedure for the three mass bins. In 

each table we indicate the bispectrum fitting /c-range (for the power spectrum the A;-range is always 0.04 — 
0.3/iMpc“^). Green plus (red minus) signs are indicating which statistics is included (excluded) from the 
fit. We see that adding the tree level bispectrum significantly improves the constraining power of the fitting 
procedure. 


4.2 Determining from the matter power spectrum 

Before we proceed with determining all the bias coefficients, the EFT c<j(i) dark matter parameter 
needs to be evaluated. This coefficient enters the EFT calculation of the halo-matter power spectrum 
equation (31) and is determined by fitting the one-loop EFTofLSS dark matter power spectrum to 
nonlinear data from iV-body simulations. As shown in [2], the one-loop EFTofLSS power spectrum 
takes the form 


k‘^ 

Feft-i-1oop = .Pii + .Pi-ioop — 2(27r)c^(-;^)-^Pii , (62) 

^NL 

and is expected to be in good agreement with data up to a scale of k O.ShMpc Beyond 
this scale, it would not be consistent to neglect the two-loop contribution [3]. By htting equation 
(62) to the power spectrum measured from the Al-body simulations (see sec. 4.1) in the fc-range 
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fcmax [h/Mpc] 

binO 

binl 

bin2 

mm 

0.22 

-0.31 

0.22 - 

0.31 

0.22 - 

0.31 

hm 

0.24 

-0.35 

0.22 - 

0.35 

0.22 - 

0.50 

hh 

0.19 

-0.32 

0.17- 

0.30 

0.20- 

0.35 

mmm 

0.14 

-0.22 

0.14- 

0.22 

0.14- 

0.22 

hmm 

0.13 

-0.22 

0.13- 

0.22 

0.13- 

0.21 

hhm 

0.13 

-0.22 

0.13- 

0.22 

0.14- 

0.23 

hhh 

0.13 

-0.21 

0.13- 

0.21 

0.11 - 

0.23 


Table 1: Table of estimated /Cmax values for the power spectra and bispectra. Power spectra values 
correspond to the range where theoretical error estimates obtained from higher perturbative order 
are estimated to be about 3%. Bispectra fcmax ranges correspond to the scale where theoretical errors 
obtained from higher perturbative order are expected to reach the same value as the variance of 
simulation data given by equation (61). 

[0.15 — 0.3]hMpc“^ one hnds: 



(2.31 ±0.02) X 


1 


^NL 

hMpc“^ 


2 


(63) 


In Fig. 1 is shown the resulting comparison of the one-loop EFTofLSS power spectrum divided by 
the nonlinear spectrum Pnl- As detailed in [10], we also place uniform 2% errorbars on the data to 
account for possible systematic errors in the simulation and in the comparison procedure between 
the simulation and the EFT. The residual oscillations of about 2% in size around the nonlinear data 
originate from bulk flow effects and can be appropriately taken into account by resumming the IR 
modes (see [13]). 


4.3 Procedure for comparing to data 

We briefly describe below the method that allows us to compare to simulation data and at the same 
time to determine the bias parameters from the same data, in particular from power spectrum and 
bispectrum measurements. We employ a total of hve observables: halo-matter (cross) and 

halo-halo (auto) power spectra as well as halo-halo-halo halo-halo-matter and 

halo-matter-matter bispectra. The EFT predictions for these observables involve all together 

seven independent bias parameters: 65 ^ 2 , ^ 5 , 3 , ^ 52 , bse and a constant contribution which we 

label Constg. These bias coefficients correspond to our choice of the basis {BoD). Note that only hve 
of these {bs,i, 65 , 2 , 5^2, bse and Conste) enter the bispectrum results at tree level, and bse parameter 
does not enter the one loop power spectra. 

Previous results on the dark matter one-loop power spectrum in the EFTofLSS are believed to be 
valid up to scales k ~ 0.3hMpc“^ [2]. Considering that dark matter helds are the building blocks of 
halos observables, we should thus not assume for halo cross or auto power spectrum that the fc-range 
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binO 

binl 

bin 2 

^ 5,1 

1.00 ± 0.01 

1.32 ±0.01 

1.89 ±0.02 

b5,2 

0.23 ±0.01 

0.52 ±0.01 

-0.51 ±0.04 


0.48 ±0.12 

0.66 ±0.13 

1.48 ±0.28 

bs2 

0.28 ± 0.01 

0.30 ±0.01 

1.88 ±0.03 

be. 

0.72 ±0.16 

0.27 ±0.17 

-0.66 ±0.37 

bse 

0.31 ±0.08 

0.76 ±0.17 

1.38 ±0.30 

Conste 

5697 ± 108 

10821 ± 169 

28856 ± 570 


Table 2: Best fit bias parameters table. After renormalization seven bias parameters remain to be 
evaluated. We use Wbody simulations to evaluate the best fit values for these parameters on hve 
observables: halo-matter cross power spectrum Phm, halo-halo auto power spectrum Phh, the two 
cross halo-matter-matter and halo-halo-matter bispectra Bhmm, Bhhm and dually the auto halo-halo- 
halo bispectra Bhhh- Parameter values are given for three mass bins (see [25, 24] for the mass ranges), 
binO being the lightest and bin2 the most massive. 


is any larger. In order to reduce the cosmic variance effect on the two point statistics, instead of using 
the halo-halo and halo-matter power spectra directly, we divide these by the dark matter 

predictions so as to obtain the ratios: = pi^h) jp{mm) ^ ^{hm) _ p{hm) jp{mm) 1 - 22 ^ 33 ] jn 

this way one is able to reduce the error in simulations due to sampling variance at low fc’s, where 
the linear approximation becomes exact. 

We complement the above by adding uniform 2% errorbars on the quantities and 
obtained from the simulations data to account for possible systematics. We also note again here that 
all the hnal bias components in eq. (31) and in eq. (32) (i.e. the functional forms multiplying the 
different bias terms) take into account the appropriate resummation of IR modes, as shown in [13] 
and displayed for the dark matter case in Fig. 1. 

We construct the halo-halo Xhh halo-matter x\.m function, which depend polynomially on six 
bias parameters enters only in the bispectra): 



E 




(^y)i 


ki<k„ 


r(k) 




(64) 


where the {xy) indices may represent both, halo-matter (hm) cross-correlation or halo-halo {hh) auto¬ 
correlation. The theoretical prediction is evaluated at the ki points of simulation measurements 
denoted as . The error on the spectrum ratios as already noted above, consists of two 

components: the leading error estimate is derived from the equation (60), to which a 2 % is added to 
account for possible unkown systematics in simulations and in the comparison between the EFTofLSS 
and simulations. Finally, fcmax,r(fc) is the largest wavenumber in the power spectrum ratio considered 
in the sum. 
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k[h/Mpc] k[h/Mpc] 

Figure 3: Halo-matter (upper panels) and halo-halo (lower panels) power spectrum ratios are shown for 
the mass binO {bs = 1.0) on the left, and for binl (bs = 1.33) on the right. The solid blue line is our 
theoretical prediction containing seven bias parameters (see table 2 ), divided by the nonlinear simulation 
data measurements. The thin gray lines correspond to the theoretical error estimates stemming from our 
neglecting two-loop corrections. The horizontal thin dashed gray lines signal the 2% error that we input to 
account for the systematics in simulations. Consistently, we see that the theory stops matching the data 
when the theoretical error becomes sizable. 


In a similar fashion, one may construct the function for the three bispectrum observables: 

1 2 


^^hxy 




B^^y{kii, k2j, ksi) - B’^^j^{kii, k2j, k^i] 


ma.'x.,B(k) 




(65) 


where x and y can represent halo h or matter m . The theoretical bispectrum predictions B^^^ are 
expressed according to equation (22) and evaluated at the ki points of the simulation measurements 
B^. Again, the errors of the bispectra ^B^^ also have two components: the hrst is given by the 
simulation error estimate in equation (43), the second is the 2% covering for possible unidentihed 
systematics. The quantity fcmax,s(A:) is the largest wavenumber allowed into the sum 

^^Incidentally, we remind the reader that the bispectrum depend only on three bias parameters 65^1, 65^2, 

and 652, which is, at tree level, equivalent to the local in time approximation (see e.g. [ 36 , 34 , 35 , 22 ] ), but this 
degeneracy is expected to break down at higher orders. On the other hand, and include two shot noise 

contributions that were not considered earlier in the literature, even in the local-in-time approximation: Const^ and 
b^s- Furthermore, we remind the reader that all former analytic techniques, when applied to halos, are harmed by the 
fact that the analytic description of dark matter is incorrect as soon as loop corrections are considered. 
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Assuming independence for both the two and three point statistics, i.e. neglecting the correlation 
between any of the functions, the total Xtot ^ of components: 


2 2,2 
Xtot - Xrhm + 


rhh 


+ xL 


+ xL 


+ xL, 


( 66 ) 


Since Xtot ^ polynomial of sixth order in the bias parameters, it has well dehned minima (assuming 
reasonable constraints on parameters). As a consequence, minimization can be achieved by applying 
any of the standard minimization tools (as found e.g. in the Mathematica package). 

This is a good point to discuss our choice of the maximum wavenumber over which, for each given 
observable, we perform the £t to data. We choose this interval with care in order to avoid as much 
as possible any issue with overfitting. First of all, we can estimate the /c-reach of each observable 
by using the estimates in sec. 2.5. The results are reported in table 1 We therefore choose the 
following fcmax for the various observables and bins over which we optimize the £t between theory 
and simulations. For the power spectra, for each bin the fc-interval is shown in Fig.s 3 and 7 below as 
delimited by two vertical lines between approximately 0.03hMpc“^ < fcfit < 0.26hMpc“^. Instead, 
for the bispectra, we choose the fcmax depending on the bin: /Cmax,B = 0.16, 0.11, O.lOh Mpc“^, for 
bin0,l, and 2 respectively. It is important to notice two aspects of these choices. First, the risk of 
overhtting by choosing too a large /cmax in the bispectrum is very limited once the available number 
of data is large enough. This is due to the complexity of the functional form of the bispectrum 
For the power spectra, where there appear three bias coefficeints that do not appear in the bispectra, 
the risk of overhtting is more realistic, and we therefore limit the fc-region of the fit to a region which 
is consistent with our expectations from the theoretical errors. Second aspect to notice is that in 
order to establish the estimates for the fc-reach in table 1, one needs to have the numerical values 
of the biases, which can be determined only after htting. We have therefore implemented the htting 
procedure iteratively, obtaining a self-consitent solution. 


4.4 Results and Plots 

In this part we present bias parameters measured from the simulated data of halo-matter, halo-halo 
power spectra as well as halo-halo-halo, halo-halo-matter and halo-matter-matter bispectra. 

In Fig. 2 we show the p-values and minimal for individual and joined best hts. We start from 
joint hts just for the two point functions, we do so for the three different halo mass ranges. We 
find that using seven free parameters gives good p-value up to the theoretically expected fc-ranges 
for both, power spectra and bispecta statistics. Adding the three bispectrum observables results in 
more constraining power and gives a p-value of 30% for binO (6yi ~ 1.0) with a maximal scale of 
fcmax,B = 0.16hMpC"h 

notice in passing that this error bar is determined by when the higher order terms are estimated to be a 3% 
fraction of the non-linear results, after subtracting for the shot noise, which dominates the signal at relatively high fc’s. 

Unfortunately, the fc-reach of bin2 is so low that we have only about 32 bispectra data. This is dangerously small 
a number, and with a dangerously large error due to the large cosmic variance at low k’s, to sufficiently determine all 
our bias coefficients with confidence. This means that there is a risk that the results for bin2 are overfitted, and they 
should be therefore taken with great caution. This issue should not apply to binO and binl, where the bispectra data 
points we use are respectively 112 and 41. 
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Figure 4: p-values corresponding to three halo mass bins, shown as a function of the maximum triangle 
side length /cmax,B for the bispectra (and fixed fc m ax p = 0.3/iMpc“^). We see the characteristic sharp 
drop in the p-value after the maximal scale fcmax,_B- As expected for higher mass bins results fail earlier, 
around fcmax,^ = 0.08 for binO and then at slightly higher scales for middle mass halos (around scales of 
km£Lx,B = 0.09) for binl and bin2). For lithest mass binO expansion preforms up to scales ~ fc ma x P = 0.15) 

Increasing the halo mass (i.e. higher bs^i) whilst keeping the p-value in the range 15 — 40%, causes 
the scale to drop to the value O.llhMpc"^ and O.lOhMpc"^ for binl and bin2 respectively. 

We also note that, having the same value of fc max r for binl and 2 is due to the discreteness of our 
simulation data points, and that a hner k binning would probably slightly increase the range fc ma x R 
of binl or/and decrease it for bin2 (as suggested by the p-value in table 2). 

In table 2 we show the best £t values of the bias parameters for the three halo mass bins. The 
parameters are obtained by means of a global £t to all five observables in Fig. 2 (last line in the 
tables) with the corresponding fcmax,B and with fcmax,p = 0.3hMpc“^. The error bars of the best £t 
parameters given in the table are obtained by estimating the 50% drop in the of p-value relative to 
the one given by the best fit parameters. 

In Fig. 3 we show the results for the halo-matter and halo-halo power spectrum ratios for mass binO 
(on the left), and binl (on the right), using the parameter values given in table 2. We show the (ratio 
of) power spectra ratios relative to the results obtained from iV-body simulation measurements. As a 
shaded region in the plot we also display the estimate for the theoretical error due to our neglecting 
two-loop effects, multiple and divided by 2 to account for our uncertainties. As always, the 2% 
error corresponding to the potential systematics in the simulations in included. Notice that in the 
cross correlation between matter and halos, the deviation of the prediction of the EFTofLSS from 
the non-linear data happens approximately when the theoretical error becomes relevant. This is 
an important conhrmation of the correctness of the EFTofLSS, whose predictions fail when they 
are expected to fail. The failure of the EFTofLSS happens at the fc-values as estimated from the 
theoretical error also for the halo-halo power spectrum, even though this is not visible in the plot, 
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because the signal is dominated by shot noise, while the theoretical error is computed relative to 
the subleading non-shot noise signal. This is why the theoretical error appears small. Similar plots 
corresponding to the higher halo mass bin2 are in Fig. 7 in Appendix B. The plots of all the bispectra 
(cross and auto): halo-matter-matter, halo-halo-matter, halo-halo-halo, for all the mass binO can be 
found in the Appendix Fig. 8 , 9, 10. 

In Fig. 4 we show the p-values corresponding to three halo mass bins, as a function of the 
maximum triangle side length fcmax,B for the bispectra (and hxed fcmax.p = 0.26hMpc“^). We see the 
characteristic sharp drop in the p-value after the maximal scale fc ma x R. Within these scales our EFT 
results £t the data well, and then fail after crossing of r scales. Higher mass (higher bias 6 yi) 
objects have typically smaller fc max r range. The sharp drop in the p-value as a function of fcmax,s is 
indication of the power of the bispectra in constraining the parameters of the theory. 

Finally we look at the results when we exclude bias parameters b^s or 652 from our theory pre¬ 
dictions. We hnd that removing the e6 operator induces a sharp drop in the p-value. This is not the 
case if one alternatively removes the operator in 5^. This suggests that the functional form of some 
of these parameters is similar, given the error bars of our data. We leave a study of the similarity 
of the contribution of the various operators to future work. We notice that once the number of data 
points is much larger than the number of parameters, it is not a big issue to have degenerate, or 
quasi degenerate, parameters, unless one is using large prefactors to force unnatural cancellations 
between the contributions, as it is not the case here. 

4.5 Configuration-space results for the cross and auto power spectra 

For completeness,we show the results in conhguration space. Correlation functions can be obtained 
by simple Fourier transform 

— = j ^ P^^{k)jo{kr) , (67) 

where jo is the spherical Bessel function. The only subtle point in this procedure is that the EFTofLSS 
in Fourier space fails at high fc’s beyond the estimated theoretical reach. When it fails, it does so with 
large deviations from the non-linear power spectrum. In order for the errors at high wavenumber not 
to erroneously affect the correlation functions at large distances, we extrapolate the power spectra 
in the EFT to high k beyond their reach, about k ~ 0.3hMpc“^ for one loop and k ~ 0.6hMpc“^ 
for two loops, by using some well behaved function. Results are insensitive to the procedure when 
the extrapolation function is reasonable. For example, we extrapolated the power spectrum by 
continuously connecting the EFT result with a rescaled linear power spectrum or by applying a 

Gaussian damping e without finding appreciable changes in the plots we show below. 

In Fig. 5 we show the dark matter cross correlation. In Fig. 6 we show the cross and auto 
correlation functions corresponding to the halo-halo and halo-matter power spectra. For both figures, 
the Icr error bars on the simulation data are given by the square root of the diagonal elements of 
the respective covariance matrices calculated in the Gaussian limit following [49]. Focussing on the 
region of the BAG peak, we can see that IR-resummation of the contribution from the displacements 
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Figure 5: Dark matter correlation function. The subscript ir means that we have performed the IR 
resummation. Error bars from simulations represent the estimated 2cr (see [49]). Error bars are correlated 
though we do not quantify their correlation. 


is crucial to reach agreement with the non-linear power spectrum (see for example [9]). Subsequent 
non-linear corrections at one loop and two loops (for dark matter) seem to provide a relatively small 
effect in that region, shifting the curve within the simulation error bars. At short distances, the 
theory curves begin to fail, as expected by the fact that the power spectrum in the EFT stops being 
correct at high wavenumbers. Given that at low wavenumbers the EFT power spectra do not present 
any apparent residual, we indeed find that at long distances, and in particular around the BAO peak, 
the results are compatible with the error bars from simulations. 

5 Conclusions 

Understanding the behavior of the Large Scale Structures (LSS) of the Universe has become a 
vital prerequisite if we wish to continue to gather information about the primordial universe. The 
EFTofLSS provides a novel technique to analytically describe the LSS at large distances, where 
perturbations are small. So far, the EFTofLSS has been compared with simulation data for statistics 
involving the dark matter field, and, as we have described in the Introduction, the success has been 
remarkable. In this paper we have first extended the theoretical description of biased tracers in the 
EFTofLSS by including the effect of baryonic physics and primordial non-Gaussianity. We have then 
compared the predictions for biased tracers of the EFTofLSS against dark-matter-only simulations 
with Gaussian initial conditions. We have focused on five different statistics: the halo-halo power 
spectrum, the halo-matter cross correlation at one-loop level, and the halo-halo-halo, halo-halo- 
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Figure 6 : Halo-matter and halo-halo correlation functions. Error bars from simulations represent the 
estimated 2a region (see [49]). Error bars are correlated though we do not quantify their correlation. The 
green curve is obtained by using only the 65 ^ bias term and linear theory, with tree-level IR-resummation, 
while the blue curve represent the EFTofLSS prediction at one-loop, with one-loop IR-resummation. 


matter, and halo-matter-matter bispectra at tree level. We have found that the predictions of the 
EFTofLSS for biased tracers at this order depend on seven bias parameters, after taking account of 
degeneracies. Four of these parameters enter in all of the hve statistics we consider. According to each 
mass bins, we have estimated the theoretical error associated to the next order terms, allowing us 
to estimate the fc-reach of each computations, and finding that lower mass bins should have a larger 
fc-reach than the higher ones. We have restricted our study to three mass bins, and found that, 
at the order at which we compute, the two-point functions from EFTofLSS match the numerical 
data to percent level up to A; ~ 0.3hMpc“^, while the bispectra match until k ~ 0.16hMpc“^, 
k ~ 0.11hMpc“^ and k ~ 0.10hMpc“^ for binO, binl and bin2 respectively. These results are 
consistent with the theoretical expectations We also shown results for the same two-point functions 
as just described, but in configuration space. Ending consistent results. 

We believe that our findings, if confirmed by other studies and comparisons, open the way to un¬ 
derstanding in an analytical way the statistics of biased tracers up to quite larger wavenumbers than 
pervasively believed, and to complement the role of simulations at smaller distances where physics 
is non-linear. Hopefully, this complementary approach will allow us to extract much cosmological 

^^Notice that our simulations have ag = 0.9, implying that the universe we compare with is quite more non-linear 
than the one we currently live in. This suggests that the fc-reach of our calculations will be slightly larger in more 
realistic cosmologies. 
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information from on-going and next generation LSS surveys. 
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Appendix 

A Q operators 

In this appendix we provide a direct connection between the operators in equation (9) and the ones 
identihed by in eq. (10): 


1 st order: 


( 68 ) 


Cg(k,i) = ^<‘'(k,«). 


C,(k,t) = [e]k, 


c«i(k,() = M]L‘', 


2 nd order: 


Cg(k,i) = [a,iW 

cg(k, i) = i®(k, t) - , 

Cg,(k,i) = rt">(t), 

Cih(k.i) = y]b(i). 

Cg(k,i) = l€®]k. 


3rd order: 


Cg(k,i) = t 
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cg(k,() = - [a^p) i 


+ [S.S;'5“’ ) . 


Cg\(k,i) = -2|i<‘>8.i'‘>^9<‘>]k(«), 


9* 


82 


cf> (k, 0 = [il=»(() + , 


-'52,2 


82 


cg,(k,i) = rt’>w, 

clf,(k,i) = -2|sWa.(s'’")<‘'^9''>]k(i). 

CS^k, i) = + 2h“S.(s"“)'‘)^9l‘>]k(i). 

Cgi(k,i) = Mr(i). 

C2,’i(k,() = |si|k’((), 

C 5 ?’i(k,i) ^^‘“’(k.t), 

Cg>^,(k,i) = |fc'l=l(() , 

Below we give the explicit form of the Cj operators introduced in the equation (10) in term of linear 
densities and kernels: 

^f\Kt) = J (69) 

^f\Kt) = j ^i(?i,?2)5§^(k-qi-q2)5^^)(qi,t)<5(^nq2,^), 

Cf^(k,t) = j ^^l{q^,q2,qs)6^o\k-q^-q^-q^)6^^\q^,t)6^^\q2,t)6^^\q^,t), 

C^Sjk,t) = 6if{k,t), 

where n-th kernels are symmetrized in the momentum variables. The unsymmetrized version of 
the kernels is given by: 


1 st order: 

2 s!i(qi) = 1 . 


( 70 ) 


2 nd order: 
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4^(qi>q2) 
45(^1’^ 2 ) 
4?i(qi>q2) 
^J!i(qi,q 2 ) 


(qi,q 2 ) - ^ 

1, 

(qrq 2 )^ _ 1 
qhl 3’ 


3rd order: 
4^(qi>q2,q3) 

4!5(qi)q2,q3) 

45(qi)q2,q3) 

4M(qi,q2,q3) 

^/2yqi,q2,q3) 

43li(qi’q2>q3) 

^J^i(qi,q2,q3) 

^jyqi,q2,q3) 

^J?i(qi,q2,q3) 

ciM(qi>q2)q3) 

4a(qi)q2,q3) 

4!^i(qi>q2>q3) 


1 f (qi-q2+qrq3) r'(2) /„ n '1 -u qi-q2(qi-q3+q2-q3) A 

2 {ql+ql+ 2 nr^ 3 ^ iq 2 ,q 3 J+ 52^32 

(F( 2 ) (q^, q 2 ) g 2 _ q^ . q^) , 

(qi+q2)-q3 

2 ^fg| 


F^^^ (qi, q2, q 3 ) + ^""^2^212"^' (qi ■ q2 - 2^(2) (q^, q2) q^) 


_ qi-(q2+q3) _^( 2 ) / \ 

2(5|+g|+2q2.q3)‘-^ lq 2 , q3J 


_2 q2q3 

q'i 


2 F( 2 ) (q^,q 2 ) + 2 a| 


qs 


oHaJk f (qrq2)^ 

93 V <ih2 


2 f’®(qi.q 2 )( 


((qi+q2)'q3)^ 

(qi+q2)^'?3 


^ +2 


£1233 

93 


(£ 1132 )^ 

q'Ul 


9qrq2qrq3q2-q3-3(qrq3)^9|-3(qrq2)^g3+g?(-3(q2-q3)^+2gig3) 

99?9l9| 

(G'(') (qi, q 2 ) - (q^, q^)) (^iag|)! _ i j 


(qi, qs, qa) - F^^^ (qi, qs, qs) 

+ 2 F( 2 ) (q^, q 2 ) (F( 2 ) (q^ + q 2 , qg) - (q^ + q^, qg)) 

(qi-q 2 )^ _ 1 

q'iq'i 3 


There are degeneracies amongst the set of operators in equation (68). Using the BoD basis 
introduced in the section 2.1 we are able to identify degenerate operators and express them as linear 
combinations of the basis element: 

2nd order operators: (71) 

(p( 2 ) _ _ 17(p(2) 

^^ 2,1 — 2^5,2 6^52,15 


3rd order operators: 

(^<(3) _ _Zfr'(3) 17^^(3) 

^^ 2,1 — 2^(5 ,2 6 ^52,1’ 
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Cl 

CAi 


21^(3) 

I 9 


''.■^7 _1_ ^(p(3) 

'( 5,2 ' /I '‘-A'? 


103^(3) 


4 '^6,2 4 ^^5,3 

7r'(3) I 9r’(3) _ 37„ , , 

2 ^(5,2 “T 2 ^(5,3 12 ^ 52^1 

_ 9(p(3) _|_ 9(p(3) _ (p(3) _ ^ 

— ^^S,3 ^,52,1 42 

7(p(3) I 7(p(3) _ U 

0 ^x 2 1 “T a^X2o fi 


's 2 , 2 ’ 


137 (p(3) I 511 (p(3) _ 3(p(3) 

16 ^(52,2 “*■ 72 ^< 53,1 4^«2 . 

'(p(3) 71(p(3) I 25(p(3) l(p(3) 

'^^21 24^2^2 'T' 12^3,1 2 ^ 2 , 2 ’ 

'(p(3) _|_ (p(3) 2(p(3) 

'^52_2 ^ 53^1 7^„2 . 


"* 2 ^ 2 ’ 


■ ' H- -^L-' ' — — 

'(52,1 4^(52,2 6 ^(53,1- 


Finally, we repeat here the expression for the bias tracer overdensity field given in terms of the BoD 
operator basis (eq. (12)) more explicitly: 


4(k, t) = cyi(t) (cg(k, t) + cg(k, t) + cg(k,t)) (72) 

+ Cy2(2) (i) (k, i) + C(5,2(3) (^) (k, i) 

+ cys (t) c(kc) + Cy3., {t) (k, f) 

+ C2 l('2)(t) C^2^^(k, t) + C(52,1(3)(^) C^2^j^(k, t) 

+ Cs 2 ^ 2 {t) cj^^k, t) + C32^2(^) cj^k, t) 

+ C53^i(t) C§yk,t) 

+ Ce 4 (f) C[^i(k, f) + ... 


where the new bias coefficients q are given in terms of the old ones from equation (10): 


C(5,l = cs,i, 

C(5,2{2) = Cy2 + ics2,l , 

C(5,2(3) = Cy2 ~ 2^2,1 + xC3,l + 2^1,1 + 2c^,i , 

C(5,3 = Cy3 ^C3,l + |Csl,l + 2c^^i , 


Cy3cs 

C52^1(2) 

5524 ( 3 ) 

£ 52,2 

Cs 2,2 

C53,l 


^S,3cs ’ 

17 

^ 17^ 103^ 145^ 97^ I 21^ 

g g ^s^,l 25 4 ; 

1 37 

^S ^,2 

3 12 

^s ^,2 4^5^,! 

^ I 511^ 25^ I ^ 17^ 

^(5^,1 ~T~ 72 ^5^,1 6 • 


(73) 


At last, we present an explicit expression for the third term of eq. (18): 

d?q 


3 Pii{k]t,t) 


3 ^f^k,-q,q)^ii(g;t,t) = 


/ (277 

C54(f) 3 Pii(A;;f,f) 


AD 


(277) 


c,i(k,-q, q)^ii(gCC) 
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+Cs: 


3 Pii{k;t,t) j (k,-q,« 

[, q)Pii(g;t,t) 
,,q)Pii(g;t,t) 


+C5,2(3)(^) 

+C5,3(t) 3 Pii(fc;t,t) y -q, 

,i(3)W 3 Pii(fc;t,t) j ^^c^jyk,-q^ 

4;^ c^jyk,-q,q)Pii(g;tp) 


+C53 i(t) 3 Pii 


+0^2 2(t) 3 Pii 


+Cs2 ,2(^) 


and for the second term in eq. (20): 

^44 q,q)] Vii(g;t,t)Pn(|k- q| 




3Pn(fc;t,t) y ^ 


\-,t,t) = 

Pii(g;t,t)Pii(|k- q|;t,t) 


+C5,i 

+C5,2(2)(t)c52_l(2)(t) 4 y 44 


)42i(k-q, q)Ai( 
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/ 44 42 (k-q,q: 

+c| 2 ( 2 )(^) 2 y 44 ^ii(p^>^)^ii(|k-q|;^,^) 

+C52,i(2)(t) 2 y 44 42i(k-q,q) ^ii(?;^,^)^ii(|k-q|;t,t) 
+C5,i(t)c5,2(2)(t) 4 J 44 ^^i^k - q> q)42(k - q, q)^ii(?; t, t)Pii(|k - q|; t, t) 
i{t)cs, 2 C 2 )it) 4 y 44 ^ q)Ai(?;^P)^ii(|k - q|;tp) 

r rJ^n /^\ /r.\ 

L(g;t,t)Pii(|k-q|;t,t) 


(75) 


B Results for one additional halo mass 


In this appendix we present the rest of the power spectra (for higher mass bin2) and bispecta plots. 
In figure 7 we show the halo-halo and halo-matter power spectrum results for this higher mass bin2, 
similar as in hgure 3. 

In the rest we present the figures for the bispectra halo-halo-halo, halo-halo-matter and halo- 
matter-matter bispectra as given in the equations (22) with the best fitted parameters in table 2. 
binO results are given in hgures 8, 9 and 10. binl results are given in hgures 11, 12 and 13. bin2 
results are given in hgures 14, 15 and 16. 
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Figure 7: Same as Fig. 3 for the higher halo mass bin2. 


The attentive reader will notice that some of the plots in e.g. Fig. ( 8 , 9, 10) noticeably deviate 
from unity. In interpreting these plots one should keep in mind the combined effect of the relatively 
large simulation error estimate due to cosmic variance (black dashed lines) and the fact that we model 
the bispectrum only up to tree level, and so the fc-reach is limited. For the sake of completeness, in 
these plots we show a number of bispectrum conhgurations even though for these we are well outside 
the fc-reach of our model, whilst for the corresponding longer scales predicted by the EFT there are 
no available simulation data. For these particular conhgurations then, a conspicuous deviation from 
unity is to be expected. In fact, as shown in the main text, the p-value of the ht of the theory to 
data is very large up to the fc-reach of the theory. 
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Figure 8: Halo-matter-matter bispectrum results for binO {bsi = 1.0). All the lines are divided by the 
simulations data. Blue line is the tree level theory prediction. Black dashed lines represent the simulation 
error bars as given in equation (61), and gray lines represent the estimate of the theory errors given by the 
equations (43). 
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Figure 9: Halo-halo-matter bispectrum results binO {bsi = 1.0). Lines are the same as in figure 8. 
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Figure 10: Halo-halo-halo bispectrum results binO {bsi = 1.0). Lines are the same as in figure 8 . 
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Figure 11: Halo-matter-matter bispectrum results binl {bsi = 1.33). Lines are the same as in figure 8. 
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Figure 12: Halo-halo-matter bispectrum results binl {bsi = 1.33). Lines are the same as in figure 8. 
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Figure 13: Halo-halo-halo bispectrum results binl (bsi = 1.33). Lines are the same as in figure 8. 
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Figure 14: Halo-matter-matter bispectrum results bin2 {bsi = 1.95). Lines are the same as in figure 8. 
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Figure 15: Halo-halo-matter bispectrum results bin2 {bsi = 1.95). Lines are the same as in figure 8. 
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Figure 16: Halo-halo-halo bispectrum results bin2 (bsi = 1.95). Lines are the same as in figure 8. 
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